library(Seurat)
library(scRepertoire)#
library(ggplot2)
library(cowplot)
library(scater)
library(scran)
library(dplyr)
library(Matrix)
#library(muscat)
library(reshape2)
library(celldex)
library(BiocParallel)
library(BiocNeighbors)
library(data.table)
library(DEsingle)
library(stringr)
library(sva)
library(readxl)
library(DESeq2)
library(DESeq)#
library(pamr)
library(ggpubr)
library(ggraph)
library(gplots)
library(pca3d)#
library(rgl)
library(scatterplot3d)
library(FactoMineR)
library(ggfortify)
library(useful)
library(tidyverse)
library(kableExtra)
library(xfun)
library(psych)
library(limma)
library(calibrate)
library(pheatmap)
library(ggraph)
library(circlize)
library(scales)
if(!require(pca3d)){install.packages("pca3d")}
# Sys.chmod("C:/Program Files/R/R-4.2.2/library",'777') #提权，手动提权
BiocManager::install("DESeq")
BiocManager::install("pca3d")
BiocManager::install("scRepertoire")
# install.packages("DESeq")
#add seurat
# setwd("~/gse162498")
setwd("C:/Users/forbing36/Desktop/单细胞生信/GSE212966")

Only_T <- readRDS(file="./data/temp/Only_T_cluster_id_test.rds") 
# 定义顺序，方便出图
levels(x = Only_T)
levels(x = Only_T) <- c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", "NKT", "CD8_CTL", "CD8_Trm", "CD8_Te")
levels(x = Only_T)
table(Idents(Only_T))
# CD4_Tn   CD4_Th CD4_Treg  CD4_Tem      NKT  CD8_CTL  CD8_Trm   CD8_Te 
# 2540       62     2811     3619      827     1109     2064     2018 
DimPlot(Only_T, reduction = "umap", label = TRUE, pt.size = 0.5) 
DimPlot(Only_T, reduction = "umap", label = FALSE, pt.size = 0.5) 
table(Idents(Only_T))

#细胞及细胞中基因与RNA数量
slotNames(Only_T)
#assay
Only_T@assays
dim(Only_T@meta.data)
View(Only_T@meta.data)

##T分群绘图####
#each type of cells
hms_cluster_id<-Only_T
hms_cluster_id<-readRDS("./data/temp/T_cluster_id_test.rds")
# CD4_Tn<-subset(Only_T, idents=c('CD4_Tn'))
# DimPlot(CD4_Tn, reduction = "umap")
# saveRDS(CD4_Tn, file="CD4_Tn.rds")


CD4_Tn<-subset(hms_cluster_id, idents=c('CD4_Tn'))
pdf("./data/output/CD4_Tn_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD4_Tn, reduction = "umap")
dev.off()
saveRDS(CD4_Tn, file="./data/temp/CD4_Tn.rds")

CD8_Trm<-subset(hms_cluster_id, idents=c('CD8_Trm'))
pdf("./data/output/CD8_Trm_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD8_Trm, reduction = "umap")
dev.off()
saveRDS(CD8_Trm, file="./data/temp/CD8_Trm.rds")

CD8_Te<-subset(hms_cluster_id, idents=c('CD8_Te'))
pdf("./data/output/CD8_Te_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD8_Te, reduction = "umap")
dev.off()
saveRDS(CD8_Te, file="./data/temp/CD8_Te.rds")

CD4_Tem<-subset(hms_cluster_id, idents=c('CD4_Tem'))
pdf("./data/output/CD4_Tem_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD4_Tem, reduction = "umap")
dev.off()
saveRDS(CD4_Tem, file="./data/temp/CD4_Tem.rds")

CD4_Treg<-subset(hms_cluster_id, idents=c('CD4_Treg'))
pdf("./data/output/CD4_Treg_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD4_Treg, reduction = "umap")
dev.off()
saveRDS(CD4_Treg, file="./data/temp/CD4_Treg.rds")

CD8_CTL<-subset(hms_cluster_id, idents=c('CD8_CTL'))
pdf("./data/output/CD8_CTL_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD8_CTL, reduction = "umap")
dev.off()
saveRDS(CD8_CTL, file="./data/temp/CD8_CTL.rds")

CD4_Th<-subset(hms_cluster_id, idents=c('CD4_Th'))
pdf("./data/output/CD4_Th_DimPlot.pdf",width = 8,height = 6)
DimPlot(CD4_Th, reduction = "umap")
dev.off()
saveRDS(CD4_Th, file="./data/temp/CD4_Th.rds")

NKT<-subset(hms_cluster_id, idents=c('NKT'))
pdf("./data/output/NKT_DimPlot.pdf",width = 8,height = 6)
DimPlot(NKT, reduction = "umap")
dev.off()
saveRDS(NKT, file="./data/temp/NKT.rds")




#DEG差异分析####
#deg in CD4_Tn
# table(Idents(Only_T))
# CD4_Tn  CD8_Trm   CD8_Te  CD4_Tem CD4_Treg  CD8_CTL   CD4_Th 
# 2540     2064     2018     3619     2811     1109       62 


# #失败🌟⚠️🪲，调整，补救差异基因merge的时候好的样本放在后面，正常的样本放在后面。
# CD4_Tn<-readRDS("./data/temp/CD4_Tn.rds")
# a<-CD4_Tn@meta.data
# write.table(a,"./data/temp/a.csv",sep=",")
# 
# table(CD4_Tn@meta.data$tech)
# CD4_Tn@meta.data$new_tech <- CD4_Tn@meta.data$tech
# levels(CD4_Tn@meta.data$new_tech)
# CD4_Tn@meta.data$new_tech <- as.factor(CD4_Tn@meta.data$new_tech)
# levels(x = CD4_Tn@meta.data$new_tech)
# levels(x = CD4_Tn@meta.data$new_tech) <- c("PDAC", "ADJ")
# levels(x = CD4_Tn@meta.data$new_tech)
# 
# b<-CD4_Tn@meta.data
# write.table(b,"./data/temp/b.csv",sep=",")
# 
# #merge失败
# CD4_Tn<-readRDS("./data/temp/CD4_Tn.rds")
# CD4_Tn_ADJ <- subset(CD4_Tn, subset=CD4_Tn@meta.data$tech=="ADJ")
# 
# CD4_Tn@meta.data$tech


class(CD4_Tn)
CD4_Tn.sec<-as.SingleCellExperiment(CD4_Tn)

table(CD4_Tn@meta.data$tech)
# 正常样本
group<-factor(c(rep(1,1214),rep(2,1326)))

rds<-readRDS('./data/temp/CD4_Tn.rds')
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
# 27579  2540
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/results_CD4_Tn")
write.table(results.classified,"./data/temp/results.classified_CD4_Tn")

#deg in Natural_Killer_T
Natural_Killer_T<-readRDS("Natural_Killer_T.rds")
b<-Natural_Killer_T@meta.data
write.table(b,"b")
class(Natural_Killer_T)
Natural_Killer_T.sec<-as.SingleCellExperiment(Natural_Killer_T)
group<-factor(c(rep(1,3139),rep(2,3890)))

rds<-readRDS('Natural_Killer_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"results_Natural_Killer_T")
write.table(results.classified,"results.classified_Natural_Killer_T")

#deg in T_Helper
T_Helper<-readRDS("T_Helper.rds")
c<-T_Helper@meta.data
write.table(c,"c")
class(T_Helper)
T_Helper.sec<-as.SingleCellExperiment(T_Helper)
group<-factor(c(rep(1,1791),rep(2,1317)))

rds<-readRDS('T_Helper.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"T_Helper")
write.table(results.classified,"results.classified_T_Helper")

#deg in Effector_Memory_CD4_T
Effector_Memory_CD4_T<-readRDS("Effector_Memory_CD4_T.rds")
d<-Effector_Memory_CD4_T@meta.data
write.table(d,"d")
class(Effector_Memory_CD4_T)
Effector_Memory_CD4_T.sec<-as.SingleCellExperiment(Effector_Memory_CD4_T)
group<-factor(c(rep(1,3226),rep(2,5101)))

rds<-readRDS('Effector_Memory_CD4_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Effector_Memory_CD4_T")
write.table(results.classified,"results.classified_Effector_Memory_CD4_T")

#deg in Resident_Memory_CD8_T
Resident_Memory_CD8_T<-readRDS("Resident_Memory_CD8_T.rds")
e<-Resident_Memory_CD8_T@meta.data
write.table(e,"e")
class(Resident_Memory_CD8_T)
Resident_Memory_CD8_T.sec<-as.SingleCellExperiment(Resident_Memory_CD8_T)
group<-factor(c(rep(1,1674),rep(2,1248)))

rds<-readRDS('Resident_Memory_CD8_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Resident_Memory_CD8_T")
write.table(results.classified,"results.classified_Resident_Memory_CD8_T")

#deg in Terminally_Exhausted_CD8_T
Terminally_Exhausted_CD8_T<-readRDS("Terminally_Exhausted_CD8_T.rds")
f<-Terminally_Exhausted_CD8_T@meta.data
write.table(f,"f")
class(Terminally_Exhausted_CD8_T)
Terminally_Exhausted_CD8_T.sec<-as.SingleCellExperiment(Terminally_Exhausted_CD8_T)
group<-factor(c(rep(1,2373),rep(2,2866)))

rds<-readRDS('Terminally_Exhausted_CD8_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Terminally_Exhausted_CD8_T")
write.table(results.classified,"results.classified_Terminally_Exhausted_CD8_T")

#deg in Regulatory_T
Regulatory_T<-readRDS("Regulatory_T.rds")
g<-Regulatory_T@meta.data
write.table(g,"g")
class(Regulatory_T)
Regulatory_T.sec<-as.SingleCellExperiment(Regulatory_T)
group<-factor(c(rep(1,3675),rep(2,2160)))

rds<-readRDS('Regulatory_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Regulatory_T")
write.table(results.classified,"results.classified_Regulatory_T")

#deg in Cytotoxic_CD8_T
Cytotoxic_CD8_T<-readRDS("Cytotoxic_CD8_T.rds")
h<-Cytotoxic_CD8_T@meta.data
write.table(h,"h")
class(Cytotoxic_CD8_T)
Cytotoxic_CD8_T.sec<-as.SingleCellExperiment(Cytotoxic_CD8_T)
group<-factor(c(rep(1,289),rep(2,731)))

rds<-readRDS('Cytotoxic_CD8_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Cytotoxic_CD8_T")
write.table(results.classified,"results.classified_Cytotoxic_CD8_T")

#deg in Effector_Memory_CD8_T
Effector_Memory_CD8_T<-readRDS("Effector_Memory_CD8_T.rds")
i<-Effector_Memory_CD8_T@meta.data
write.table(i,"i")
class(Effector_Memory_CD8_T)
Effector_Memory_CD8_T.sec<-as.SingleCellExperiment(Effector_Memory_CD8_T)
group<-factor(c(rep(1,29),rep(2,966)))

rds<-readRDS('Effector_Memory_CD8_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Effector_Memory_CD8_T")
write.table(results.classified,"results.classified_Effector_Memory_CD8_T")

#deg in Pre_Exhausted_CD8_T
Pre_Exhausted_CD8_T<-readRDS("Pre_Exhausted_CD8_T.rds")
j<-Pre_Exhausted_CD8_T@meta.data
write.table(j,"j")
class(Pre_Exhausted_CD8_T)
Pre_Exhausted_CD8_T.sec<-as.SingleCellExperiment(Pre_Exhausted_CD8_T)
group<-factor(c(rep(1,238),rep(2,223)))

rds<-readRDS('Pre_Exhausted_CD8_T.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"Pre_Exhausted_CD8_T")
write.table(results.classified,"results.classified_Pre_Exhausted_CD8_T")

#volcano_plot####
setwd("~/gse162498/deg")
res <- read.csv("deg_Cytotoxic_CD8_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Effector_Memory_CD4_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Effector_Memory_CD8_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,8),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Natural_Killer_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_CD4_Tn", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Pre_Exhausted_CD8_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_T_Helper", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Resident_Memory_CD8_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Terminally_Exhausted_CD8_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-2,2),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Regulatory_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-7,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)
